LOCAL DISCONTINUOUS GALERKIN METHODS FOR 
ONE-DIMENSIONAL SECOND ORDER FULLY NONLINEAR 
ELLIPTIC AND PARABOLIC EQUATIONS 

XIAOBING FENG* AND THOMAS LEWISt 

Abstract. This paper is concerned with developing accurate and efficient discontinuous Galerkin 
methods for fully nonlinear second order elliptic and parabolic partial differential equations (PDEs) 
in the case of one spatial dimension. The primary goal of the paper to develop a general framework 
for constructing high order local discontinuous Galerkin (LDG) methods for approximating viscosity 
solutions of these fully nonlinear PDEs which are merely continuous functions by definition. In 
order to capture discontinuities of the first order derivative Ux of the solution u, two independent 
functions qi and q2 are introduced to approximate one-sided derivatives of u. Similarly, to capture 
the discontinuities of the second order derivative u^x, four independent functions pi, p2, ps, and p4 
are used to approximate one-sided derivatives of qi and q2. The proposed LDG framework, which is 
based on a nonstandard mixed formulation of the underlying PDE, embeds a given fully nonlinear 
problem into a mostly linear system of equations where the given nonlinear differential operator must 
be replaced by a numerical operator which allows multiple value inputs of the first and second order 
derivatives Ux and Uxx- An easy to verify criterion for constructing "good" numerical operators is also 
proposed. It consists of a consistency and a generalized monotonicity. To ensure such a generalized 
monotonicity, the crux of the construction is to introduce the numerical moment in the numerical 
operator, which plays a critical role in the proposed LDG framework. The generalized monotonicity 
gives the LDG methods the ability to select the viscosity solution among all possible solutions. 
The proposed framework extends a companion finite difference framework developed by the authors 
in |9] and allows for the approximation of fully nonlinear PDEs using high order polynomials and 
non-uniform meshes. Numerical experiment results are also presented to demonstrate the accuracy, 
efficiency and utility of the proposed LDG methods. 
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1. Introduction. This is the third paper in a series [5| [TU] which is devoted 
to developing finite difference (FD) and discontinuous Galerkin (DG) methods for 
approximating viscosity solutions of the following general one-dimensional fully non- 
linear second order elliptic and parabolic equations: 

(1.1) F {uxx,Ux,u,x) ^ 0, X e ^ := {a,b) CR, 
and 

(1.2) ut+F{u,.„u,,u,x,t)^0, ix,t)enT:^nx{0,T), 

which are complemented by appropriate boundary and initial conditions. 

Fully nonlinear PDEs, which are nonlinear in the highest order derivatives of the 
solution functions in the equations, arise in many applications such as antenna design, 
astrophysics, differential geometry, fluid mechanics, image processing, meteorology, 
mesh generation, optimal control, optimal mass transport, etc (cf. [51 section 5]), 
and, as a result, the solution of each of these application problems critically depends 
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on the solution of their underlying fully nonlinear PDEs. In particular, it calls for 
efficient and reliable numerical methods for computing the viscosity solutions of these 
fully nonlinear PDEs. Currently, the availability of such numerical methods are very 
limited (cf. 8J). 

The goal of this paper is to design and implement a class of local discontinuous 



Galerkin (LDG) methods for the fully nonlinear equations (1.1 1 and (1.2 1. The more 



involved high dimensional generalizations of the LDG methods of this paper will be 
reported in [TT]. 

Because of the full nonlinearity, integration by parts, which is the necessary tool 
for constructing any DG method, cannot be performed on equation (1.1 1. The first 
key idea of this paper is to introduce the auxiliary variables p := Uxx and q :— Ux and 
rewrite the original fully nonlinear PDE in the following nonstandard mixed form: 



(1.3) 
(1.4) 
(1.5) 



F{p,q,u,x) = 0, 
q-Ux ^0, 
P-Qx =0. 



Unfortunately, since Ux and Uxx may not exist for a viscosity solution u e C°(f2), 
the the above mixed form may not make sense. To overcome this difficulty, our second 
key idea is to replace q = Ux hy two possible values of Ux, namely, its left and right 
limits, and p = 9^ by two possible values for each possible q. Thus, we have the 
auxiliary variables qi,q2 : ^ R and pi,p2,P3,P4 ■ 17 —> R such that 



(1.6) 

(1.7) 

(1.8) 

(1.9) 

(1.10) 

(1.11) 



qi{x) - Uxix^ 
q2{x) - Ux{x+ 
pi{x) - qixix^ 
P2{x) ~ qix{x'^ 
Paix) - q2xix^ 
P4{x) ~ q2x{x^ 



= 0, 
= 0, 
= 0, 
-0, 
-0, 
= 0. 



We remark that (1.6) paired with the equation (1.8) or (1.9), and (1.7 1 paired with 



equation (1.10 1 or (1.11 1, can each be regarded as a "one-sided" Poisson problem in 
u (in a mixed form) with source terms pi, P2, Ps, P4, respectively. 



To incorporate the multiple values of Ux and Uxx, equation (1.3) must be modified 



because F is only defined for single value functions p and q. To this end, we need the 



third key idea of this paper, that is, to replace (1.3) by 



(1.12) 



Fipi,P2,P3,P4.,qi,q2,u,x) = 0, 



where F, which is called a numerical operator, should be some well-chosen approxi- 
mation to F. 

Natural questions now arise regarding to the choice of F. For example, what are 
criterions for F and how to construct such a numerical operator? These are two im- 
mediate questions which must be addressed. To do so, we need the fourth key idea of 
this paper, which is to borrow and adapt the notion of the numerical operators from 
our previous work [9] where a general finite difference framework has been developed 
for fully nonlinear second order PDEs. In summary, the criterions for F include con- 
sistency and g-monotonicity (generalized monotonicity), for which precise definitions 
can be found in section [2j^It should be pointed out that in order to construct the de- 
sired numerical operator F, a fundamental idea used in [5] is to introduce the concept 
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of the numerical moment, which can be regarded as a direct numerical realization for 
the moment term in the vanishing moment methodology introduced in |12| (also see 
section 4], [13]). 



Finally, we need to design a DG discretization for the mixed system (1.6)-(1.12| 
to complete the construction of our LDG method. This then calls for the fifth key idea 
of this paper, which is to use different numerical fluxes in the formulations of LDG 
methods for the four "one-sided" Poisson problems in their mixed forms described 
by (1.6) - We remark that, to the best of our knowledge, this is one of a 



few scenarios in numerical PDEs where the flexibility and superiority (over other 
numerical methodologies) of the DG methodology makes a vital difference. 

This paper consists of four additional sections. In section[2]we collect some prelim- 
inaries including the definition of viscosity solutions, the definitions of the consistency 
and g-monotonicity of numerical operators, and the concept of the numerical moment. 
In section[3]we give the detailed formulation of LDG methods for fully nonlinear ellip- 
tic equation following the outline described above. In section|4]we consider both 
explicit and implicit in time fully discrete LDG methods for fully nonlinear parabolic 



equation ( 1.2 1. The explicit four stage classical Ronge-Kutta method and the implicit 
trapezoidal method combined with the spatial LDG methods will be explicitly con- 
structed. In section [5] we present a number of numerical experiments for the proposed 



LDG methods for the fully nonlinear elliptic equation (1.1) and their fully discrete 



counterparts for the parabolic equation (1.2). These numerical experiments not only 



verify the accuracy of the proposed LDG methods but also demonstrate the efficiency 
and utility of these methods. 

2. Preliminaries. Standard space notations are adopted in this paper. For 
example, B{n), USC{^1) and LSC{il.) denote, respectively, the spaces of bounded, 
upper semi-continuous, and lower semicontinuous functions on Q. For any v G B{V,), 
we define 

w*(x) := lim sup 11(1/) and i;*(a;) := liminf u(y). 

y — y 

Then, v* G USC{Q) and e LSC{Q), and they are called the upper and lower 
semicontinuous envelopes of v, respectively. If v is continuous, there obviously holds 
V = V* — v^,. 

Let F : S"^^"^ xR'^xRxfi— >Rbca bounded function, where 5'*^'^ denotes the 
set oi d X d symmetric real matrices. Then, the general second order fully nonlinear 
PDE takes the form 

(2.1) F{D'^u,Wu,u,x) ^0 inU. 

Note that here we have used the convention of writing the boundary condition as a 
discontinuity of the PDE (cf. 2, p.274]). 

The following two definitions can be found in [Zl [31 [2] . 

Definition 2.1. Equation (2.1) is said to be elliptic if for all (q. A, a;) e R'* x 



R X r2, there holds 

(2.2) F{A,q,X,x) < F{B,q,X,x) yA,B e S"^""^, A> B, 

where A > B means that A — B is a nonnegative definite matrix. We note that when 

OF 
dA 



F is differentiable, the ellipticity also can be defined by requiring that the matrix 



is negative semi-definite (cf. 7^ p. 441]). 

Definition 2.2. A function u e i?(il) is called a viscosity subsolution (resp. 



super solution) of (2.1 1 if, for all € C (fi), if u* — (p (resp. — p) has a local 
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maximum (resp. minimum) at xq G 17, then we have 

F4D^ip{xo),Vipixo),u*{xo),xo) < 
(resp. F*{D'^lp(xq),'V(p{xo),u^{xo),xo) > 0). The function u is said to be a vis- 



cosity solution of (2.11 if it is simultaneously a viscosity suhsolution and a viscosity 
supersolution of ( |2.1[ ). 

We note that if F and u are continuous, then the upper and lower * indices can 



be removed in Definition 2.2 The definition of elhpticity imphes that the differential 



operator F must be non-increasing in its first argument in order to be elliptic. It 



turns out that elhpticity provides a sufficient condition for equation (2.1 1 to fulfill a 
maximum principle (cf. [TJ [3]). From the above definition it is clear that viscosity 
solutions in general do not satisfy the underlying PDEs in a tangible sense, and the 
concept of viscosity solutions is nonvariational. Such a solution is not defined through 
integration by parts against arbitrary test functions; hence, it does not satisfy an 
integral identity. This nonvariational nature of viscosity solutions is the main obstacle 
that prevents direct construction of Galerkin-type methods, which require variational 
formulations to start. 

The following definitions are adapted from in the case d = 1. 

Definition 2.3. 

(i) A function F : — )■ R is called a numerical operator. 

(ii) A numerical operator F is said to be consistent ( with the differential operator 
F) if F satisfies 

(2.3) liminf i^(pi,p2,P3,P4, 9i, 92, Ai, ^i) > g. A, ^), 

Pfc — >-p,fc = l,2,3,4 

(2.4) limsup F(pi,p2,P3,P4,gi,'?2, Ai,^i) < F*(p,(3', A,$), 

where F^, and F* denote respectively the lower and the upper semi-continuous 
envelopes of F. 

(iii) A numerical operator F is said to be g-monotonc if F{pi,p2,P3,Pi, qi, q2, A, ^) 
is monotone increasing in pi and p^ and monotone decreasing in p2 andp^, 
that is, F(t,;,;,t,9i,'?2,A,C)- 

We remark that the above consistency and g-monotonicity play a critical role in 
the finite difference framework established in |9^. They also play an equally critical 
role in the LDG framework of this paper. In practice, the consistency is easy to fulfill 
and to verify, but the g-monotonicity is not. In order to ensure the g-monotonicity, 
one key idea of §) is to introduce the numerical moment to help. The following is an 
example of a so-called Lax-Friedrichs-like numerical operator adapted from [S] : 

(2.5) F(pi,p2,P3,P4,gi,g2, A,0 := F{ — - — , — - — , A,0 

-i-a(pi -P2-P3 +P4), 



where a G R is an undetermined positive constant and the last term in (2.5 1 is 
called the numerical moment. It is trivial to verify that F is consistent with F. By 
choosing a correctly, we can also ensure g-monotonicity. For example, suppose F is 
differentiable and there exists a positive constant M such that 

dF 



(2.6) M > 
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Then, it is trivial to check that the Lax-Friedrichs-hke numerical operator is g- 
monotone provided that a > M. 

We conclude this section with a few remarks about the above definitions. 

Remark 2.1. (a) By the definition of the ellipticity for F, the monotonicity 
constraints on F with respect to P2 and p^ in the definition of g-monotonicity are 
natural. 

(h) By choosing the numerical moment correctly, the numerical operator F then 
behaves like a uniformly elliptic operator, even if the PDE operator F is a degener- 
ate elliptic operator. The consistency assumption then guarantees that the numerical 
operator is still a reasonable approximation for the PDE operator. 

(c) Sometimes it may not be feasible to globally bound ; however, it is suffi- 
cient to chose a value for a such that the g-monotonicity property is preserved locally 
over each iteration of the nonlinear solver for a given initial guess. 

(d) The role of the numerical moment as well as the interpretation of the numer- 
ical moment will be further discussed in Section\5~l 



3. Formulation of LDG methods for elliptic problems. We first consider 



the elliptic problem (1.1 1 with boundary conditions 



(3.1) u(a) = Ua and u(b) = Uf, 

for two given constants Ua and Uf,. 

Let {xj}'^-^Q C be a mesh for il such that xq = a and xj ~ b. Define Ij ~ 
{xj-i,Xj) and hj = Xj — xj-i for all j — 1,2,..., J, ho — ft-j+i — 0, and h — 
max.i<j<j hj . Let Th denote the collection of the intervals which form a 

partition of the domain il. We also introduce the broken iJ^-space and broken C"- 
space 

h\th) n ^'(^)' ^"c^o ■■= n ^"(^)' 

and the broken L^-inner product 

(w,w)7-ft:=^^ / vwdx VvjW £ L'^{rt). 

For a fixed integer r > 0, we define the standard DG finite element space C 
H\rh) C L\Th) as 

:= n 

len 

where Vril) denotes the set of all polynomials on / with degree not exceeding r. We 
also introduce the following standard jump notation: 

[vh{xj)] := Vh{xj) - Vh{x+) for j = 1, 2, • • • , J - 1. 



We now are ready to formulate our LDG discretizations for equations (1.6 )-( 1.12 ) 



First, for (fully) nonlinear equation ( 1.12 1 we simply approximate it by its broken L 
projection into V^, namely, 

(3.2) ao{uh,qih,q2h,Pih,P2h,P3h,Pih;(l>oh) = V^oh e V'^, 
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where 



an {uh ,qih, q2h , Pih , P2h , Psh , Pih ; '/'o/i ) = [F{pih, P2h , Psh ,P4:h,qih, q2h ,Uh,-), '/'o/i 



n 



Next, we discretize the four /meor equations (1.8)-(1.11|. Notice that for given 
{pijUi, (|L6| and Of, (fTel) and ([L9), (|1.7| aEd ([rTo|, and (ILT]) and 



sources 



(1.11) are four (different) Poisson equations for u. Thus, we can use the mixed up- 



winding LDG formulation for the Laplacian operator to discretize these equations. 
The only difference in the four equations will be how we choose our upwinding nu- 
merical fluxes for Uh, qih and q2h- To realize the above strategy, we first define the 
element-wise LDG formulation, and we then define the whole domain LDG formula- 
tion afterward. 

3.1. Element-v^rise LDG formulation. Suppose that values for Uh{a~), Uh{b'^), 
qih{0'~), and qih{h^) for « 1,2 are given. We postpone explaining how these val- 
ues are chosen until section 3.2 Our LDG discretization of equations ( |1.6[ )-( [m| ) is 
defined as follows: for all (t>ih <E , 



(3.3) 



qih (pih dx+ Uh {(j)ih)x dx = u/,.(xj) (jy^hixj ) - Uh{xj^i) (j)ih{xj_^) 



fori -1,2, J -1,2, 



and for all ^kh G , 



, J, and 



- ifi = l, 
+ ifi = 2; 



(3.4) / pkh ipkh dx+ q-f^f^ {ipkh)x dx = q-^hi^j) ^khix, ) - qi,h{Xj-i) ^kh{xj_-^) 



for k = 1,2,3,4, and 



1 if /c is odd, 

2 if /c is even, 



- if fc = 1, 



+ if A: = 2. 

Notice that the nodal values determined by a and j3 follow directly from equations 



(1.6)-(1.11) 



3.2. Boundary numerical fluxes. To complete the construction, we must 
specify how the boundary numerical flux values for Uh^ qih, and q2h are determined 
in the above formulation. Due to the inherent jumps of piecewise constant functions, 
which corresponds to the case r = 0, we shall consider the two cases r > 1 and r — 
separately. 

When r > 1, we have freedom to control how the functions u/j, qih, and q2h 
approach the boundary. Thus, we can assume continuity across the boundary for 



Uh, qih, and q2h- Considering the boundary conditions given by (3.1), the continuity 
requirement naturally leads to 



(3.5) 



Uhia''^)=Ua, Uh{b'^) = Ub. 
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On the other hand, since no boundary data for qih or q2h is given, any choice 
of the boundary numerical fluxes for them is a guess (unless one aheady knows the 
exact solution u). Here we choose 

(3.6) Qihia^) ^ qih{a'^), qih{b~^) = qih{b~), 1 = 1,2. 

It is important to note that both qih{a^) and qih{b~) {i = 1,2) are treated as un- 
knowns in the above LDG formulation. The choice ( |3.6[ ) is equivalent to requiring 
that qih is continuous at the boundary nodes x = a and x = b. 

We now consider the case r = 0. To define the boundary numerical fluxes, we 
first examine the consequences of the interior flux choices represented by the above 
LDG formulation. Suppose 7^ is a uniform mesh and denote the midpoint of Ij by 
Xj for all Ij e Th- Define Uj :— u'' (xj). Then, it follows from (3.3) and (3.4) that 

(3.7) qih (xj) 

(3.8) q2h (xj) 

(3.9) pih [Xj) - - — o^Uj-1, 



h 








h -'^-^1' 








qihixj) - qih{xj-i) 




2 - 2C/,- 


I + 


h 

qih{xj+i) - q\h{xj) 




-1 - 2C/, 


+ Uj+i 


h 

q2h{X]) - q2h{xj-i) 






+ Uj+1 


h 

q2h{X] + l) - q2h{xj) 






+ Uj+2 


h 









(3.10) (%) = = ■■= S',U„ 

, . . ^ q2h{xj) - q2h{xj-i) Uj-i - 2Uj + Uj+i 

(3.11) [Xj] = 1 = .= d^Uj, 

(3.12) Pih [Xj] - - p .- O^Uj+i, 

for j = 3, 4, . . . , J — 2. Thus, in order to define boundary values for u/j, qih, and 
q2h, we need to define ghost values C/-i, Uq, Uj+i, and Uj+2 that are equivalent to 
extending the solution u to the outside of the domain fl. Below we describe a natural 
way to do such an extension that is consistent with the interpretation of the auxiliary 
variables. 

From the Dirichlet boundary data for u, a natural choice is that Uq = Ua and 
[/j+i = Uh- This is equivalent to assuming 

(3.13) u^{a-) = Ua, u^{h+) = Ub. 

In other words, we extend the boundary data for u away from the boundary over 
an interval of length h. Due to the inherent discontinuities of the piecewise constant 
functions, u^{a^) and u^{b~) are treated as unknowns. Otherwise, the boundary data 
would be extended into the interior of the domain over an interval of length h. 

From ( |3.7[ ) and ( |3.8[ ) we see that q2h (xj) = qih (%+i) and qih (xj) = q2h (xj-i) 
in the interior of the domain. Extending this relationship to the boundary yields 

(3.14) q2h{a^) = <7i/i(a^), qih{b'^) = q2h{b^), 

where both qih{a^) and q2h{b ) are treated as unknowns in the above LDG formula- 
tion. 

Finally, we need to define values for qih{a^) and q2h{b'^)- Using 12J as a guide, 
we are led to choosing 

(3.15) qih{a^) = qih{a^), q2h{b'^) = q2h{b~). 
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We note that this is consistent with discrctizing the auxihary boundary conditions 

(3.16) (gih).(a) - (<Z2/0.(fo) =0. 

In order words, we require that qih and q2h are constant across the boundary. Using 
ghost values, the above requirements are equivalent to imposing the constraints 

U^l^2Ua-Ui, Uj+2=2ut-Uj. 

Remark 3.1. From the imposed boundary conditions we can see that the rela- 
tionship p2h — P3h has been extended to the boundaries. Thus, using the ghost values 
defined above and substituting the equations (3.7)-(3.12) into (3.2|, we successfully 
recover the convergent finite difference method defined in J^l for the grid function U . 
Thus, for r = 0, the convergence of the proposed LDG method is obtained. Heuris- 
tically, using higher order elements should increase the rate and/ or accuracy of con- 
vergence. 

3.3. Whole domain LDG formulation. Using the above element-wise LDG 
formulation (3.3 1 and (3.4 1, and substituting the boundary numerical flux values from 
section [3^ we get the following whole domain LDG discretization of ( |1.6| )-( [lTT| ): 

(3.17) [qih, (f>ih)r^ + ai{uh, (t)ih) = fi{(l>ih), ^ff^ih e V'^, « = 1, 2, 

(3.18) {Pjh^i^jh)^^ + bj{qih,q2h;ipjh) = 0, V-0jh e F'', j = 1,2,3,4, 

where 

.7-1 



J-1 

a2{v, if) = [v, (p^)r^ + (1 - Kr) w(a+) <^(a+) - ^ v{x+) [(p{xj)] , 

3 = 1 

J-1 

biivi,V2;(p) = {vi,(p^)r^ + vi{a+) Lp{a+) - vi{b')(p{b^) - ^vi{x~)[(p{xj)], 

i=i 
J-1 

b4{vi,V2;ip) = {v2,'Px)n + W2(a"^) <^(a+) - V2{b~)(p{b~) - ^V2{x'^)[(p{xj)], 

i=i 

&2(wi,W2;</?) = {vi,(px)rh +wi(a+)(^(a+) - (1 - Hr) V2ib^) ip{b^) 

J-1 

- Kr Vl{b'')ip{b^) -^Vi{x+)[ip{Xj)], 

b3{vi,V2;ip) = {v2,(px)n + (1 ~ i^r) Wi(a+) (^(a+) + U2(a+) (/?(«+) 

J-1 

- V2{b~) (p{b~) - V2{xj)[ip{x/' 



^J7J ' 



and 



fl{cj)) = KrUb(j){b )-Ua(t>ia+), 
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for 

(3.19) 



if r = 0, 

1 otherwise. 



In summary, our nonstandard LDG methods for the fully nonlinear Dirichlet 
problem ((lT|) and pi) are defined as seeking (uh,qih,q2h,Pih,P2h,P3.h,Pih) e {V^) 



such that (3.2), (3.17), and (3.181 hold. 



We conclude the section with a few remarks. 

Remark 3.2. (a) Looking backwards, (3.171 and (3.181 provide a proper inter- 
pretation for each of qih and pjh for i = 1, 2 and j = 1, 2, 3, 4, for a given function Uh- 
Each qih defines a discrete derivative for Uh and each pjh defines a discrete second 
order derivative for Uh- The functions qih and q2h should be very close to each other 
if Ux exists. Similarly, the functions pih, P2h, Psh, o,nd p^h should be very close to 
each other if u^x exists. However, their discrepancies are expected to be large if or 
Uxx, respectively, do not exist. The auxiliary functions qih defined by (3.17) and the 



auxiliary functions Pjh defined by (3.18) can be regarded as high order extensions of 
their lower order finite difference counterparts defined in 19]. 



(b) It is easy to check that the linear equations defined by (3.17|-(3.18) are linearly 
independent. 



(c) Notice that (3.2 1, (3.17), and ( 3.18^ form a nonlinear system of equations, 
with the nonlinearity only appearing in oq. Thus, a nonlinear solver is necessary in 
implementing the above scheme. In section^ an iterative method is used with initial 
guess given by the linear interpolant of the boundary data. Since a good initial guess is 
essential for most nonlinear solvers to converge, another possibility is to first linearize 
the nonlinear operator and solve the resulting linear system first. However, we show 
in our numerical tests that the simple initial guess works well in many cases. We 
suspect that the g-monotonicity of F enlarges the domain of "good" initial guesses 
over which the iterative method converges. 

4. Formulation of fully discrete LDG methods for parabolic problems. 

The goal of this section is to extend the LDG methods for elliptic problems to solving 
the initial-boundary value problem (1.2) using the method of lines. Let the initial 
condition be given by 



(4.1) 



u{0,x) = uq{x), 



and the boundary conditions be given by 



(4.2) 



u{t,a) = Ua{t), u{t,b) — Ub{t), V<e(0,r] 



We shall consider both the implicit trapezoidal rule and the (explicit) fourth order 
classical Runge-Kutta method (i.e., RK4) for the time-discretization. In practice, 
the time-discretization scheme should be chosen to match the order of the spacial 
discretization. Thus, when using a piecewise constant element, a sufficient choice for 
the time discretization would be forward or backward Euler. However, when using 
higher order elements, a higher order scheme such as RK4 should be chosen. 



We first formulate the semi-discrete in space discretization for equation (1.2) 



which is a straightforward adaptation of the one described in section [3j Let (p £ V 



Replacing the PDE operator with a numerical operator in (1.2), and using the LDG 
framework of section [3j we obtain the semi-discrete equation 

(4.3) {uht,(t>h)-j-^ ^ -{f {pih,P2h,P3h,P4h,qih;q2h,Uh,t,-) 
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where, given Uh at time t, corresponding values for qih and pjh can be found using 
the methodology below. 

We now describe a full discretization procedure for (1.2 1 by applying an ODE 
solver to the semi-discrete equations given in (4.3). For a fixed integer M > 0, 
let At ^ jj and tk := k At for k e (0, M]. Notationally, ul{x) e V'^ will be an 
approximation for u(t„, x) for n = 0, 1, . . . , M . We define the initial value to be 
the -projection of Uq, namely, 



(4.4) 



Next, we introduce several "one-sided" discrete differential operators, which will 
be used to define explicit time-stepping schemes and to define the auxiliary variables 
at time i = in implicit time-stepping schemes. We first define two "one-sided" first 
order discrete derivatives QiV, Q2V € for a given function «(•, tk) G H^{Th) by 



(4.5) (Q^, 



Kr Ub{tk) (f>l{b ) - Uaitk)4>lia'^) + (l - Kr)v{b ) (6 ) 
,7-1 



(4.6) {Q2V,(j}2)j-^ = Ubitk)<j)2ib )-K^Wa(ife)02(a+)-(l-Kr)l'(a+)(/'2(a+) 



,7-1 



The above definitions are inspired by (3.17). The super- index k on Qf indicates that 
the definitions are i-dependent because of the boundary terms. 

We also define four discrete "one-sided" second order discrete derivatives ViV, V2V, 
V^v,V^v e at time tk by 



(4.7) {V^v,^i) 



Th 



(4.8) (7'^,^4 



(4.9) (7'^,V^2) 



Q^«(6-)^i(6-)-Qjt;(a+)Vi(a+) 
,7-1 

- (2^7', ^ix)^^ + E 2^(^7) [Mx,)], 
7-1 

(1 - Kr) Q^vib-) ^2{b-) + nr 2^(6- ) 7/'2 ) 

-2^(a+)^2(a+)-(2i«,^2.)^^ 

7-1 

+ ^2^(4) [V^2(x,)], VV;2e^^ 



(4.10) {Vlv,ij3)r, = Q2^{b-)Mb-) - il-nr)Qlvia+)Ma^) 

- Kr Q2v{a^) Ma'^) - {Q2V, V'Sx), 



7-1 
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where is defined by (3.19). The above four definitions are motivated by (3.18). 
Lastly, to simphfy the presentation, we introduce the operator notation 

(4.11) Fk[v] ■.= F{V^v,V^v,V^v,r^v,Q'iv,Q2^,v,tk,x). 

Using the new notation, the semi-discrete equation can be rewritten compactly as 

(4.12) {uht{tk,-),(l)h)r, ^ -{Fk[uh{tk,-)] ,<l>h).j-,^ V0,, fc = 1,2,--- ,M. 

4.1. The fourth order classical Runge-Kutta method. A straightforward 
application of the fourth order classical Runge-Kutta (RK4) method to (4.12) yields 



where 





'l>h] 










(6, 
















n = 










'f>h] 











Notice that in the above explicit time-stepping scheme, the function is de- 
fined as an L^-projection of the source data based on u^~^ ■ However, the boundary 
conditions are not enforced in the definition for u^. To take care of the boundary 
conditions, we choose to enforce them weakly, which requires the introduction of a 
modified L^-projection. Specifically, for any v £ L^(ri), we recall that the standard 
i'^-projection V^v G of v is defined by 



(4.13) 



For any v e C°{Th), wc introduce a modified L^-projection P^; : L'^{n)nC°{Th) ^ V' 
at time tk G (0, T] by 

(4.14) (P^, 0,Or, + {Kv{a)M<^) + Kvib)Mb) 



Clearly, the boundary conditions (4.2| are weakly enforced in (4.14) via a penalty 
technique, an idea which dates back to Nitsche [TT . 

Using the above discrete differential operators Q^, Vj, Pf^, and P/i for i = 1,2 and 
j — 1,2,3,4, and using the notation given in (4.11), our fully-discrete RK4 method 
for the initial-boundary value problem (1.2), (4.2), and (4.1) is defined as follows: for 
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n = 1,2,. 


.,M, 






(4.15) 






2^2""' +24"" 


(4.16) 


(^n — 1 

?i ~ 




2^1 J' 


(4.17) 


(^n — 1 




(4.18) 


?3 — 




2^2 J> 


(4.19) 


(^n — 1 






(4.20) 


7/° - 






We remark that it is 


easy to verify that the value ^4 ^ 



actually already take into 
account the boundary conditions at time t„ because the evaluation calls Q" and Q2 ' 
which in turn have the boundary conditions built-in. Thus, in the above formulation, 
the boundary condition enforcement can actually be successfully relaxed by replacing 
PJJ with P/i in (4.15). However, for other explicit methods a weak boundary condition 
enforcement method such as the above modified L^-projection is necessary, especially 
if the boundary conditions are not consistent with the initial condition. For example, 
in the forward Euler method, defined by 



(4.21) 



AtK, 



,n-l 



we can see that the approximation at time t„ relies upon the modified L^-projection 
in order to see the Dirichlet boundary condition at the current time. 



4.2. The trapezoidal method. Applying the trapezoidal rule to (4.121, we 
obtain 



and n — 1, 2, . . . ,M. Thus, using the trapezoidal rule to discretize (4.12), and using 
the implicit equalities 



and 



Plh — 'lPlh^ P2h — '2P2hi Psh — 'sPsh^ Pih — 'iPih 



for 71 = 1, 2, ... , M, the fully discrete trapezoidal LDG method for approximating so- 
lutions to ([r2|, (|4l]), and (|42]) is defined by seeking « , g^/i ' 92/i > Pi/i . -^2/1 > -P3/1 ' PJh ) ^ 
{V'y such that 



At: 



hlifOh 



n 



At 



, 71—1 TP r ll 



(4.22) 

(4.23) {q'^,,cp,h)r, +adK.<t>^h)=9^{t'',c^^h) 

(4.24) (pj^,^',/.)^^ +b, {qlh^q^H-.i^jh) = 0, V7A,„ € = 1,2,3,4, 



n 
1,2, 
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where < = ¥hUo, g^^ = Q%1 for z = 1, 2, p^"^ = V^u^ for j = 1, 2, 3, 4, and 

.7-1 

J-1 



j-i 

&iK,z;2";<^) = K,(^.K +<(a+)(/.(a+)-<r)(^r)-5]<(x-)[^(a;,)], 

.7-1 

j-i 

63K, W2 ; = t/'^)?^. + (1 - Kr) <(a+) </?(a+) + Kr {a+) (^(a+) 

.7-1 

gi (r , 0) = 77fc (r ) 0(r ) - [r )<P{a+), 



Again, is defined by (3.19). Notice that the above fully discrete formulation 
amounts to approximating a nonhomogeneous fully nonlinear elliptic equation at each 
time step using the LDG method defined in section |3] 

5. Numerical experiments. In this section, we present a series of numerical 
tests to demonstrate the utility of the proposed LDG methods for fully nonlinear 



PDEs of the types (1.1) and (1.2). In all of our tests we shall use uniform spatial 
meshes as well as uniform temporal meshes for the time-dependent problems. To 
solve the resulting nonlinear algebraic systems, we use the Matlab built-in nonlinear 
solver fsolve for the job. For the elliptic problems we choose the initial guess as the 
linear intcrpolant of the boundary data. For parabolic problems, we let ul = Vhuo, 
and then define all auxiliary variables by 9°^ — Q^u^^ for i = 1, 2 and = 'Pj'^h 
j = 1,2,3,4. Also, the initial guess for fsolve at the nth time step will be chosen as 
the computed solution at the previous time step when using implicit methods. The 
role of the numerical moment will be further explored in section [5. 3| 

For our numerical tests, errors will be measured in the L°° norm and the 
norm, where the errors are measured at the current time step for the time-dependent 
problems. For both elliptic and parabolic test problems where the error is dominated 
by the spatial discretization errors, it appears that the spatial errors are of order 
C'(/i'"+^). Thus, the schemes appear to exhibit an optimal rate of convergence in both 



5.1. Elliptic test problems. We first present the results for four test problems 
of type (1.1 ). Both Monge- Ampere and Bellman types of equations will be tested. 
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Test with N = 32. r= 0. and 0= 10 Test with N= 32, r= 1,and □ = 10 
1 1 1 1 6 I , , , 




• u_N * u_N 

-0,1 I 1 1 1 1 ^ H -0.1 I ^ ^ ^ 1 - 

Q 0.2 0.4 Q.6 Q.8 1 0.2 O.'^ 0.6 0.8 1 



r 


Norm 


h = 1/4 


h = 


1/8 


h = 


1/16 


h = 


1/32 






Error 


Error 


Order 


Error 


Order 


Error 


Order 





L' 


7.1e-02 


3.5e-02 


1.02 


1.4e-02 


1.30 


7.5e-03 


0.92 




L°° 


1.3e-01 


8.7e-02 


0.57 


5.3e-02 


0.73 


2.9e-02 


0.87 


1 


L' 


1.6e-02 


5.0e-03 


1.67 


1.3e-03 


1.90 


3.4e-04 


1.95 




L°° 


2.2e-02 


6.3e-03 


1.84 


1.6e-03 


2.00 


3.9e-04 


2.00 


2 


L' 


3.1e-13 


3.0e-13 


0.03 


3.0e-13 


-0.01 


3.1e-13 


-0.01 




L°° 


7.4e-13 


6.1e-13 


0.28 


6.7e-13 


-0.14 


7.1e-13 


-0.08 



Fig. 5.1. Test 1 with a = 10. 



Test 1. Consider the elliptic Monge- Ampere problem 

-u^^ + 1 = 0, 0<x<l, 
u(0) = 0, "(1) = ^. 
It is easy to check that this problem has exactly two classical solutions: 

u^{x) — -x^, (x) — — x"^ +x, 

where m+ is convex and u~ is concave. Note that u+ is the unique viscosity solution 
which we want our numerical schemes to converge to. In section |5.3| we shall give 
some insights about the selectiveness of our schemes. 

We approximate the given problem for various degree elements (r — 0,1,2) to 
see how the approximation converges with respect to h. Note, when r = 0,1, the 
solution is not in the DG space V'^. The numerical results are shown in Figure 5.1 



We observe that the approximations using r — 2 are almost exact for each mesh size. 
This is expected since m+ G when r = 2. 

Test 2. Consider the problem 

-ul^+Uxa: + S{xf - S{x) =0, -1<X<1, 

u{-l) = - sin(l) - 8 cos(0.5) + 9, u(l) = sin(l) - 8 cos(0.5) + 9, 



where 

S{x 



J ^ cos(x2) - 4x2sin(a:|a;|) + 2cos(f ) + 2, x ^ 0, 
~ [-4a;2sin(a;|.T|) + 2cos(f) + 2, x = 0. 
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T^stwltli N = 64. r= 0, and 0= 6 Testwitli N = 64, r= 1,and □ = 






Norm 


h = 1/2 


h = 


1/4 




h = 


1/8 


h = 


1/16 


h = 


1/32 






Error 


Error 


Order 


Err 




Order 


Error 


Order 


Error 


Order 







1.8 


Q.Oc-Ol 


1.04 


4.3c 


01 


1.05 


2.1C-01 


1.04 


l.Oc-01 


1.02 




L°° 


2.5 


1.2 


0.99 


6.2c 


01 


1.00 


3.1C-01 


1.00 


1.6C-01 


1.00 


1 




2.9e-01 


6.3C-02 


2.20 


1.9c 


02 


1.73 


7.0C-03 


1.44 


2.8C-03 


1.34 




L°° 


2.9e-01 


6.4G-02 


2.17 


2.0c 


02 


1.66 


7.6G-03 


1.42 


3.1C-03 


1.30 


2 




5.7G-03 


8.2C-04 


2.80 


1.3c 


04 


2.66 


3.2e-05 


2.03 


9.1e-06 


1.81 






2.0C-02 


3.1G-03 


2.70 


4.2c 


04 


2.87 


5.5e-05 


2.94 


8.0e-06 


2.77 


3 




8.8C-04 


7.7C-05 


3.51 


3.0c 


06 


4.68 


1.4e-07 


4.42 


l.Oe-08 


3.76 






2.1C-03 


1.4C-04 


3.90 


8.6c 


06 


4.01 


5.6C-07 


3.94 


9.5C-08 


2.57 



Fig. 5.2. Test 2 with a = 6. 



This problem has the exact solution u{x) — sin {x\x\) — 8 cos (|) + + 8. Note that 
this problem is not monotone decreasing in Uxx, and the exact solution is not twice 
differentiable at x = 0. However, the derivative of F with respect to u^x is uniformly 
bounded. The numerical results are shown in Figure |5.2[ As expected, we can see 
from the plot that the error appears largest around the point x = 0, and both the 
accuracy and order of convergence improve as the order of the element increases. For 
finer meshes, we see the rates of convergence begin to deteriorate. Theoretically, we 
expect less than optimal rates of convergence due to the low regularity of the solution. 

Test 3. Consider the stationary Hamilton- Jacobi-Bellman problem with finite 
dimensional control set 



min {—9uxx + Ux — u + Six)} = 0, — 1 < a; < 1, 
e(x)e{i,2} 

u(-l) = -l, «(1) = 1, 



where 



S{x) 



-Ux"^ - A\x\^ + x\x\^ , x<Q 
24a;2 - 4|a;|3 a; > 0. 



This problem has the exact solution u{x) = x\x\^ corresponding to 9* {x) = 1 for a; < 
and 9* (x) = 2 for x > 0. Approximating the problem using various order elements, 
we have the following results recorded in Figure |5.3[ Again, due to the low regularity 
of the solution, we expect the rates of convergence to be affected. However, we still 
see increased accuracy for high order elements. 

Test 4. Consider the stationary Hamilton- Jacobi-Bellman problem with infinite 
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T^stwltli N = 64. r= 0, and 0=4 



Tesi:witriN = 64, r=1,and d = 4 









Norm 


h = 1/2 




h = 


1/4 




h = 


1/8 


h = 


1/16 


h = 


1/32 






Error 


Err 




Order 


Err 




Order 


Error 


Order 


Error 


Order 







S.Oc-Ol 


3.4c 


01 


0.53 


1.3c 


01 


1.37 


6.1C-02 


1.12 


4.4C-02 


0.49 




L°° 


8.5C-01 


5.7c 


01 


0.57 


3.5c 


01 


0.72 


2.0C-01 


0.79 


l.lc-01 


0.86 


1 




1.4C-01 


4.3c 


02 


1.72 


9.7c 


03 


2.16 


2.7C-03 


1.85 


7.3C-04 


1.87 






3.6C-01 


1.1c 


01 


1.74 


3.0e 


02 


1.87 


7.8e-03 


1.94 


2.0C-03 


1.97 


2 




2.8C-02 


3.2c 


03 


3.10 


4.0c 


04 


3.00 


5.1e-05 


2.99 


6.4e-06 


2.99 






4.0C-02 


5.5c 


03 


2.84 


7.3c 


04 


2.93 


9.3e-05 


2.97 


1.2e-05 


2.98 


3 




9.4C-03 


1.3c 


03 


2.91 


1.6c 


04 


3.01 


1.9e-05 


3.01 


2.4e-06 


3.01 






l.lc-02 


1.5c 


03 


2.91 


1.9c 


04 


3.02 


2.3C-05 


3.01 


2.9C-06 


3.01 



Fig. 5.3. Test 3 with a = A. 



dimensional control set 



inf 

o<e(x)<i 



Uxx + ^ X fJ-x -\ — u + S{x) 



0, 



1.2 <x < 4, 



u(1.2) = 1.441nl.2, m(4) = 161n4, 



where 



S{x) 



4 ln(2;)2 + 12 ln(a;) + 9 - Sa;^ ln(a;)2 - 4x^ \n{x) 



4x3 [2 ln(2;) + 1] 

This problem has the exact solution u{x) — x^lnx corresponding to the control func- 
tion 0* {x) = 2x^[2\n{xyhi] • Approximating the problem using various order elements, 



we obtain the results recorded in Figure |5.4| We can see that the approximations ap- 
pear to reach a maximal level of accuracy of about 5.0e-7 in both L^- and L°°-norm, 
which is consistent with the results in Test 1 corresponding to r = 2. 

5.2. Parabolic test problems. We now implement the proposed fully discrete 
RK4 and trapezoidal LDG methods for approximating fully nonlinear parabolic equa- 
tions of the form (1.2). While the above formulation makes no attempt to formally 



quantify a CFL condition for the RK4 method, for the tests we assume a CFL con- 
straint of the form At = Kth'^, and note that the constant Kt appears to decrease 
as the order of the element increases. Below we implement both the RK4 method 
and the trapezoidal method for each test problem. However, we make no attempt to 
classify and compare the efficiency of the two methods. Instead, we focus on testing 
and demonstrating the usability of both fully discrete schemes and their promising 
potentials. For explicit scheme tests, we record the parameter Kt, and for implicit 
scheme tests, we record the time step size At. Note that the row 0* in the figures 
corresponding to the RK4 method refers to elements with ?■ = that use the standard 
projection operator in (4.15). 
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T$stwithN = 64, r=0, and d=4 



Test with N = 64. r= l.ard 0=4 







Norm 


h = 2.S/4 


h = 


2.8/8 


h = 


2.8/16 


h = 


2.8/32 


h = 


2.S/64 






Error 


Error 


Order 


Error 


Order 


Error 


Order 


Error 


Order 







5.2 


3.3 


0.67 


1.5 


1.08 


e.lc-Ol 


1.34 


2.6C-01 


1.25 






5.7 


3.6 


0.65 


1.9 


0.91 


1.1 


0.84 


5.7C-01 


0.91 


1 




2.6C-01 


8.6C-02 


1.60 


2.6C-02 


1.72 


7.4e-03 


1.83 


2.0C-03 


1.90 






3.3C-01 


l.lc-01 


1.56 


3.5C-02 


1.67 


l.le-02 


1.71 


3.4C-03 


1.65 


2 




2.6G-03 


3.9e-04 


2.77 


6.6e-05 


2.55 


1.4e-05 


2.26 


3.2C-06 


2.09 






7.3e-03 


l.Oc-03 


2.85 


1.4e-04 


2.91 


1.9e-05 


2.81 


4.1C-06 


2.25 


3 




6.4G-05 


4.2e-06 


3.93 


3.1e-07 


3.75 


1.2e-07 


1.35 


1.2C-07 


0.08 






2.7C-04 


2.1C-05 


3.72 


1.4C-06 


3.84 


8.7e-07 


0.72 


8.8C-07 


-0-01 










Fig. 5.4 


Test 4. 


with OL = 


4. 









Test 5. Let VL ~ (0, 1), Ua{t) = t^,Ub = | + t^, and uo{x) — \x^. We consider 



the problem \\.2\ , ( |4.1[ ), and ( |4.2[ ) with 



2- +^ 



1. 



It is easy to verify that this problem has a unique classical solution u{x, t) = 0.5 4- 
+ 1. Notice that the PDE has a product nonlinearity in the second order derivative. 



The numerical results for RK4 are presented in Figure 5.5 and the results for the 
trapezoidal method are shown in Figure |5.6[ As expected, RK4 recovers the exact 
solution when r = 2. 

Test 6. Let Q = (0,2), Uait) = 1, = e^t'+i), and uo{x) = e^. We consider the 



problem ( |1.2| , ( |4.1[ ), and ( |4.2[ ) with 
and 

S{x,t) 



: \n[uxx + l) + S{x,t), 



(a;-(t + l)ln((i+l)2e(*+i)" + 



=(t+i)a 



It is easy to verify that this problem has a unique classical solution u{x,t) = 
Notice that this problem is nonlinear in both u^x and Ux- Furthermore, the exact 
solution u cannot be factored into the form u{x, t) = G{t) Y{x) for some functions 
G and Y. Results for RK4 are given in Figure |5.7[ and results for the trapezoidal 



method are shown in Figure 5.8 We note that RK4 was unstable without using the 
very restrictive values for Kt recorded in Figure |5.7[ However, for RK4, we observe 
optimal rates of convergence in the spacial variable while the rates for the trapezoidal 
method appear to be limited by the lower rate of convergence for the time-stepping 
scheme. 
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T^stwltli N = 32. r= 0. and 0= 2 at time T= 1 Test with N= 32, r= 2, and □= 2 at time T= 1 




0.2 0.4 O.S 0.8 1 0.2 0.4 0.6 0.8 1 



r 


Norm 


h = 1/4 


h = 1/8 


h = 1/16 


h = 1/32 


Error 


Error Order 


Error Order 


Error Order 







9.9e-02 
2.2e-01 


6.4e-02 0.64 
1.4e-01 0.64 


3.6e-02 0.81 
8.0e-02 0.81 


1.9e-02 0.92 
4.3e-02 0.89 


1 




5.7e-03 
8.0e-03 


1.5e-03 1.98 
2.0e-03 1.99 


3.7e-04 1.99 
5.1e-04 1.99 


9.2e-05 1.99 
1.3e-04 1.99 


2 




2.4e-08 
3.6e-08 


2.4e-08 0.00 
3.7e-08 -0.03 


2.4e-08 0.00 
3.7e-08 -0.01 


2.4e-08 -0.00 
3.7e-08 -0.01 



Fig. 5.5. Test 5: Computed solutions at T = 1 using Kt = 0.001, a = 2. 



Test 7. Let D, = (0, 27r), Ua{t) = 0, Uf, = 0, and Mo(a;) — sin(x). We consider the 
problem Ol), (IXll), and with 



Fiuj.^, Uj., u, t,x) = — min \ Ag u^x ~ cix, t) cos(i) sinfa;) — sin(t) sin(a;) >, 
e{t,x)e{i,2} L J 

where Ai — 1,A2 = i, and 



c(x, t) 



1, if < t < I and 0<x<7TOT^<t<TT and tt < x <27t, 
h, otherwise. 



It is easy to check that this problem has a unique classical solution u{x, t) = cos(i) sin(a;). 
Notice that this problem has a finite dimensional control parameter set, and the op- 
timal control is given by 



e*{t,x) = 



1, if c{x, t) = 1, 

2, ifc(a:,t) = i. 



The numerical results are reported in Figure [5?9| for RK4 and in Figure [5.10| for the 
trapezoidal method. 

Test 8. Let n ^ (0,3), Ua{t) ^ e"*, = 8e~*, and uq(x) ^ \x - Ip. We 
consider the problem (1.2), (|4.1|), and (4.2) with 



F{uxx,Ux,u,t,x) = ~ inf 

-i<e{t,x)< 



^[\x - l\uxx + 0ux} - \x - (|x-l|-3) e-*. 



It is easy to verify that the problem has the exact solution u(t,x) — \x — Ipe *. 
Notice that the above operator is not elliptic for a; = 1. Also, this problem has a 
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T^stwith N = 32. r= 0. and 0= 2 at time T= 1 Test with N= 32, r= 2, and □= 2 at time T= 1 




0.2 0.4 O.B 0.8 1 0.2 0.4 0.0 0.8 1 



r 


Norm 


h = 1/4 


h = 1/8 


h = 1/16 


h = 1/32 


Error 


Error Order 


Error Order 


Error Order 







5.9e-02 
1.9e-01 


3.4e-02 0.78 
l.le-01 0.79 


1.9e-02 0.84 
5.9e-02 0.90 


l.Oe-02 0.91 
3.0e-02 0.95 


1 




5.0e-03 
l.lc-02 


1.4e-03 1.86 
2.8e-03 2.00 


3.6e-04 1.94 
7.1e-04 2.00 


9.1e-05 1.97 
1.8e-04 2.00 


2 


^00 


l.lG-07 
1.5e-07 


l.le-07 0.00 
1.6e-07 -0.04 


l.le-07 0.00 
1.6e-07 -0.01 


l.le-07 0.00 
1.6e-07 -0.00 



Fig. 5.6. Test 5: Computed solution at T = 1 using At = 0.001 and a = 2. 



Tesiiwitli N = 32, r= 0. and D=4attimeT= 0.5 Testwiiili N= 32. r= 3, and = 4 ai: time T= 0.5 




r 


Norm 


h = 2/4 


h = 2/8 


h = 2/16 


h = 2/32 


Error 


Error Order 


Error Order 


Error Order 





L°° 


6.9e+00 
l.Oe+01 


5.7C+00 0.28 
7.9C+00 0.40 


4.1e-f00 0.47 
5.6e-h00 0.50 


2.6e+00 0.65 
3.7e+00 0.60 


0* 


L°° 


2.8e-)-00 
7.8e-h00 


1.5C+00 0.89 
5.0e-h00 0.65 


8.6e-01 0.82 
2.9e-h00 0.76 


5.1e-01 0.76 
1.6e-h00 0.85 


1 


L°° 


4.8e-01 
8.4e-01 


1.2e-01 2.05 
2.3e-01 1.87 


3.0e-02 1.93 
5.8e-02 1.99 


8.2e-03 1.88 
1.5e-02 1.92 


2 


L°° 


3.7e-02 
5.9C-02 


7.8e-03 2.23 
l.Oe-02 2.53 


1.9e-03 2.03 
2.2e-03 2.19 


4.8e-04 2.00 
5.2e-04 2.11 


3 


L°° 


l.le-03 
2.5e-03 


7.4e-05 3.95 
1.6e-04 3.93 


4.7e-06 3.99 
l.le-05 3.86 


3.01e-07 3.96 
7.66e-07 3.84 



Fig. 5.7. Test 6: Computed solutions at time T = 0.5 using Kt = 0.005,0.001,0.0005,0.0001 
for r = 0, 1, 2, 3, respectively, and a = 4. Left figure uses the standard projection operator. 
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Test With N = 32, r= 0. and D=4attlmeT= 0.5 Test with N= 32. r= 3, and o = ^m\m^J= [i.5 




r 


Norm 


h ^ 2/4 


h = 


2/8 


h = 


2/16 


h ^ 


2/32 






Error 


Error 


Order 


Error 


Order 


Error 


Order 





L'' 


2.8e+00 


1.5e+00 


0.89 


8.6C-01 


0.82 


5.1e-01 


0.77 




L°° 


7.8e+00 


5.0e+00 


0.65 


2.9e+00 


0.76 


1.6e+00 


0.85 


1 


L2 


3.8e-01 


1.3e-01 


1.62 


4.5e-02 


1.49 


1.5e-02 


1.60 




L°° 


1.3e+00 


4.0e-01 


1.74 


l.lc-01 


1.85 


3.0e-02 


1.91 


2 




2.7e-02 


6.7e-03 


2.04 


1.9e-03 


1.82 


5.3e-04 


1.85 






l.Oe-01 


1.5e-02 


2.77 


2.1e-03 


2.80 


5.4e-04 


1.99 


3 




l.le-03 


7.2e-05 


3.89 


1.3e-05 


2.46 


1.2e-05 


0.09 






5.5e-03 


4.0C-04 


3.76 


2.7e-05 


3.92 


1.3e-05 


1.05 



Fig. 5.8. Test 6: Computed solutions at time T = 0.5 using At = 0.005 and a = A. 



bang-bang type control with the optimal control given by 

e*{t,x) -- 




We can see from the results for the RK4 method in Figure 5.11 and the results for the 
trapezoidal method in Figure |5.12| that the spatial rates of convergence appear to be 
limited by the low regularity of the solution, while the accuracy appears to increase 
with respect to the element degree. 

5.3. The role of the numerical moment. In this section, we focus on un- 
derstanding the role of the numerical moment. We first give an interpretation of the 
numerical moment, and then we explore the utility of the numerical moment. The 
role of the numerical moment can heuristically be understood as follows when the 
numerical moment is rewritten in the form 

,2 fPih - P2h - P3h + Pih 
^ "I 

Letting r = 0, we see that pi-i^^^p^i^^p^'^+p-"^ ig an approximation to 

we can see that the numerical moment acts as a centered difference approximation for 
A m multiplied by a factor that tends to zero with rate 0{h^). Thus, at the PDF level, 
we are in essence approximating the fully nonlinear second order elliptic operator 
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Test With N= 32.r= 0,and o= 2 attime T= 3.1 T$stwitn N = 32, r= 3, and □= 2attlmeT= 3.1 




r 


Norm 


h = tt/2 


h = tt/4 


h = tt/8 


h = 7r/16 


Error 


Error Order 


Error Order 


Error Order 







1.5e+00 
9.6e-01 


1.3e+00 0.24 
7.7e-01 0.31 


7.1e-01 0.82 
4.9e-01 0.65 


3.3e-01 1.12 
2.9e-01 0.78 


0* 




1.2e+00 
9.2e-01 


6.9e-01 0.78 
5.8e-01 0.67 


3.0e-01 1.19 
3.2e-01 0.85 


1.4e-01 1.16 
1.8e-01 0.83 


1 




2.7e-01 
1.9e-01 


7.6e-02 1.81 
6.6e-02 1.52 


2.0e-02 1.94 
1.7e-02 1.97 


5.0e-03 1.99 
4.2e-03 2.00 


2 




7.2e-02 
6.8e-02 


1.8e-02 1.98 
1.8e-02 1.93 


4.5e-03 2.02 
4.3e-03 2.04 


l.lc-03 2.01 
l.le-03 2.03 


3 




8.3e-03 
7.6e-03 


5.7e-04 3.87 
5.1e-04 3.92 


3.6e-05 3.97 
3.2e-05 4.00 


2.2e-06 4.02 
1.9e-06 4.04 



Fig. 5.9. Test 7: Computed solutions at time T = 3.10 using Kt = 0.05,0.005,0.001,0.0005 for 
r = 0, 1, 2, 3, respectively, and a = 2. Left plot uses the standard projection operator. 



by the quasilinear fourth order operator Fp, where 

{."^XXXX 7 '^XX 7 '^X ; ^5 P ^xxxx ^ {'^XX : ^X 1 ^1 ■ 

In the Umit as p — ?> 0, we heuristically expect that the solution of the fourth or- 
der problem converges to the unique viscosity solution of the second order prob- 
lem. Using a convergent family of fourth order quasilinear PDEs to approximate a 
fully nonlinear second order PDE has previously been considered for PDEs such as 
the Monge-Ampere equation, the prescribed Gauss curvature equation, the infinity- 
Laplace equation, and linear second order equations of non-divergence form. The 
method is known as the vanishing moment method. We refer the reader to [131 [5] for 
a detailed exposition. 

We now show that the proposed schemes using a numerical moment can eliminate 
the numerical artifacts that arise as consequences from using a standard discretiza- 
tion, and in certain cases when the numerical artifacts are not fully eliminated, the 
algebraic system has enough structure to design solvers that emphasize the mono- 
tonicity in ^^"^^^ and are consistent in searching for solutions at which the nonlinear 
PDE problem is elliptic. We again consider the Monge-Ampere type problem from 



Test 1 in section 5.1 The given problem has two classical PDE solutions it+ and u 



However, there exists infinitely many functions that satisfy the PDE and bound- 



ary conditions almost everywhere, as seen by fi in (5.1). These almost everywhere 



solutions will correspond to what we call numerical artifacts in that algebraic solu- 
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Test with N= 32.r= 0, and o= 2 attime T= 3,1 Testwitn N = 32, r= 3, and □= 2attlmeT= 3.1 




r 


Norm 


h = 7r/2 


h = 7r/4 


h = 7r/8 


h = 7r/16 


Error 


Error Order 


Error Order 


Error Order 







1.2e+00 
9.2e-01 


6.9G-01 0.78 
5.8e-01 0.67 


3.0e-01 1.19 
3.2e-01 0.85 


1.4e-01 1.16 
1.8e-01 0.83 


1 




2.3e-01 
2.1e-01 


7.5e-02 1.63 
6.6e-02 1.66 


2.0e-02 1.89 
1.7e-02 1.95 


7.9e-03 1.36 
5.8e-03 1.56 


2 


L°° 


6.9e-02 
6.5C-02 


1.8C-02 1.93 
1.8C-02 1.87 


4.7e-03 1.95 
4.5e-03 2.00 


2.5e-03 0.90 
1.9e-03 1.21 


3 




8.1e-03 
7.5e-03 


5.9e-04 3.79 
5.2e-04 3.87 


l.le-04 2.42 
7.6e-05 2.77 


l.le-04 0.03 
7.3e-05 0.06 



Fig. 5.10. Test 1: Computed solutions at time T = 3.10 using At = 0.031 and a = 2. 



tions for a given discretization may correspond to these functions. It is well known 
that using a standard discretization scheme for the Monge- Ampere problem can yield 
multiple solutions, many of which are numerical artifacts that do not correspond to 
PDE solutions 18]. For example, let ^ g H^{Q, 1) \ C2(0, 1) be defined by 



(5.1) 




for X < 0.5, 
i, for a; > 0.5. 



Furthermore, suppose Xj = 0.5 for some j — 2, 3,..., J — 1. Then, when using 
a standard central difference discretization, fi corresponds to a numerical solution, 
yielding a numerical artifact. 

We now consider our discretization that uses a numerical moment. When a = 0, 
we have no numerical moment. As a result, we have numerical artifacts as in the 
standard central difference discretization case. Suppose r — 0. Then, for a 7^ 0, 



inspection of (3.9)-(3.12) yields the fact that p2h cannot jump from a value of 1 to 
a value of —1 when going across Xj = 0.5. Thus, the numerical moment penalizes 
discontinuities in pji^, j — 1, 2, 3, 4, and, as a result, the numerical moment eliminates 
numerical artifacts such as /i. Similarly, for r — 1, we can see that fj, does not 
correspond to a numerical solution. However, in this case, the algebraic system does 
have a small residual that may trap solvers such as fsolve. Thus, for r — and 
r = 1, the numerical moment penalizes differences in pjh, j — 1,2,3,4, that arise 
from discontinuities in Uh, qih-, and q2h- Hence, it eliminates numerical artifacts such 
as /i. 

We now consider r > 2, in which case /i G V^. Furthermore, since fi G , we 
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T^stwltli N = 32. r= 0. and 0= 2 at time T= 1 Test with N= 32, r= 3, and □= 2 at time T= 1 




r 


Norm 


h = 3/4 


h = 3/8 


h = 3/16 


h = 3/32 


Error 


Error Order 


Error Order 


Error Order 





L'' 
L°° 


2.1e+00 
2.1e+00 


1.4e+00 0.51 
1.5e+00 0.44 


5.3e-01 1.44 
8.5e-01 0.86 


2.9e-01 0.88 
4.9e-01 0.81 


0* 


L2 

L°° 


8.6e-01 
1.8e+00 


8.8e-01 -0.02 
1.3e+00 0.53 


5.9e-01 0.58 
7.3e-01 0.80 


2.9e-01 1.00 
3.9e-01 0.92 


1 


L°° 


2.0e-01 
3.9e-01 


7.3e-02 1.45 
l.Oe-01 1.92 


1.3e-02 2.44 
2.7e-02 1.95 


3.2e-03 2.06 
6.8e-03 1.98 


2 


L°° 


l.le-01 
l.le-01 


2.0e-02 2.49 
2.0e-02 2.42 


5.2e-03 1.98 
5.2e-03 1.99 


1.3e-03 2.02 
1.3e-03 2.00 


3 


L2 


3.0e-02 
2.9e-02 


6.8e-03 2.13 
6.9e-03 2.08 


1.7e-03 2.00 
1.7e-03 2.02 


4.2e-04 2.03 
4.2e-04 2.02 



Fig. 5.11. Test 8: Computed solutions at time T = I using Kt = 0.05,0.005,0.001,0.0005 for 
r = 0, 1, 2, 3, respectively, and a = 2. 



will end up with Uh = /i, qih = l2h, and pih = pjh for j — 2,3,4 being a numeric 
solution, where 



X-' 1 



, , , ^ : for x < 0.5, , , , I 1, for X < 0.5, 

qihix) ^< % f ^ n ^ ^"""^ P^^^'^> = 1 1 f ^nK 
I— a;+|, for X > 0.5, 1—1, for x > 0.5. 

Thus, by the equalities of pji^ j = 1, 2, 3, 4, the numerical moment is always zero and 
we do have numerical artifacts. These equalities are a consequence of the continuity 
of u/i, qihj and q2h- With the extra degrees of freedom for r > 2, we allow to be 
embedded into our approximation space , thus creating possible solutions with a 
zero valued numerical moment. The numerical moment acts as a penalty term for 
differences in pjh, j = 1, 2, 3, 4, which are a consequence of discontinuities in qih and 
q2h that naturally arise for nontrivial functions when r = or r = 1. 

Even with the possible presence of numerical artifacts for the above discretization 
when r > 2, the numerical moment can be exploited at the solver level. We now 
present a splitting algorithm for solving the resulting nonlinear algebraic system that 
uses the numerical moment to strongly emphasize the fact that the viscosity solution 
of the PDE should preserve the monotonicity required by the definition of cUipticity. 



Algorithm 5.1. 

(1) Pick an initial guess for Uh- 
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T^stwltli N = 32. r= 0. and 0= 2 at time T= 1 Test with N= 32, r= 3, and □= 2 at time T= 1 




r 


Norm 


h = 3/4 


h = 


3/8 


h = 


3/16 


h = 


3/32 






Error 


Error 


Order 


Error 


Order 


Error 


Order 







8.6e-01 


8.8C-01 


-0.02 


5.9C-01 


0.58 


2.9e-01 


1.00 






1.8e+00 


1.3e+00 


0.53 


7.3e-01 


0.80 


3.9e-01 


0.92 


1 




2.0C-01 


7.3e-02 


1.45 


1.3e-02 


2.44 


3.2e-03 


2.05 




L°° 


3.9C-01 


l.Oc-01 


1.92 


2.7e-02 


1.95 


6.8e-03 


1.98 


2 




l.le-01 


2.0e-02 


2.49 


5.2e-03 


1.98 


1.3e-03 


2.02 




L°° 


l.le-01 


2.0e-02 


2.42 


5.2e-03 


1.99 


1.3e-03 


2.00 


3 




3.0e-02 


6.8e-03 


2.13 


1.7e-03 


2.00 


4.2e-04 


2.03 






2.9e-02 


6.9e-03 


2.08 


1.7e-03 


2.02 


4.2e-04 


2.02 



Fig. 5.12. Test 8: Computed solutions at time T = 1 using At = 0.001 and a = 2. 



(2) Form initial guesses for qih, q2h, Pih, P2h, P3h, andpnu using equations (3.171 



and (3.18) 



(3) Solve equation 

(4) Solve equation ( |3.17 ) fori = 1,2 and the equation formed by averaging (3.181 
for j = 2, 3 for Uh, qih, and q2h 



(5) Solve equation ( |3.18[ ) for j = 1,2,3,4 for pih, P2h, Psh, andpih- 

(6) Repeat Steps 3-5 until the change in P^''+P3h sufficiently small. 

For the next numerical tests, we will show that using Algorithm |5.1| with a suf- 
ficiently large coefficient for the numerical moment destabilizes numerical artifacts 
such as n and steers the approximation towards the viscosity solution of the PDE. 
Let u{x) — I . Then, u is the secant line formed by the boundary data for the given 
boundary value problem. We now approximate the solution of the Monge- Ampere 
type problem from Test 1 in section [S^T] by using 100 iterations of Algorithm [O] fol- 
lowed by using fsolve on the full system to solve the global discretization given by 



(3.2), (3.17), and (3.18). We take the initial guess to be 



(0) 3 1_ 



where, for r = 0, u° is first projected into V'^. From Figure 5.13 we see that the 
numerical moment drives the solution towards the viscosity solution u+ when r = 
and a is positive. From Figure |5.14[ we see that the numerical moment also drives 
the solution towards the viscosity solution u+ when r = 2 and a is positive, despite 
the presence of numerical artifacts. From Figure [5. 15[ we see that the moment drives 



the solution towards the viscosity solution of F{uxx, u^, u, x) 



1 given by u 
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Fig. 5.13. Left: a = 40, /i = 1/40, and r = 0. Right: a = 0, h = 1/40, and r = 



Test With 20. 2, and 0-20 Test with N - 20. 2, and □ - 




Fig. 5.14. Left: a = 20, ft = 1/20, and r = 2. Right: a = 0, h = 1/20, and r = 2 



for r = and r — 2 when a is chosen to be negative. In each figure, the middle 
graph corresponds to fi. Clearly, we recover the numerical artifact corresponding to /i 
when a = 0. Thus, the numerical moment plays an essential role in either eliminating 
numerical artifacts at the discretization level or handling numerical artifacts at the 
solver level. 

We make one final note about using the iterative solver given by Algorithm |5.1[ 
Note that using fsolve to solve the full system with the initial guess given by uj^ 
resulted in either not finding a root for many tests (r — 1) or converging to a numerical 
artifact with a discontinuous second order derivative at another node in the mesh 
(r = 2). In order to use fsolve for the given test problem, the initial guess should 
either be restricted to the class of functions where P2h f^nd p^h preserve the ellipticity 
of the nonlinear operator, the initial guess should be preconditioned by first using 



fsolve with r = 0, 1, or the initial guess should be preconditioned using Algorithm 5.1 



When using ?' = and a non-ellipticity-preserving initial guess, solving the full system 
of equations with fsolve st ill h as the potential to converge to u~ even for a > 0. 



The strength of Algorithm 5.1 is that it strongly enforces the requirement that F is 
monotone decreasing in p2 and p^ over each iteration. Thus, a sufficiently large value 
for a drives the approximation towards the class of ellipticity-preserving functions if 
the algorithm converges. 
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